Extreme events in the dispersions of two neighboring particles under the influence of 

fluid turbulence 
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We present a numerical study of two-particle dispersion from point-sources in 3D incompressible 
Homogeneous and Isotropic turbulence, at Reynolds number Re ~ 300. Tracer particles are emitted 
in bunches from localized sources smaller than the Kolmogorov scale. We report the first quantita- 
tive evidence, supported by an unprecedented statistics, of the deviations of relative dispersion from 
Richardson's picture. Deviations are due to extreme events of pairs separating much faster than 
average, and of pairs remaining close for long times. The two classes of events are the fingerprint of 
complete different physics, the former being dominated by inertial subrange and large-scale fluctua- 
tions, while the latter by the dissipation subrange. A comparison of relative separation in surrogate 
white-in-time velocity field, with correct viscous-, inertial- and integral-scale properties allows us to 
assess the importance of temporal correlations along tracer trajectories. 
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The relative separation of pairs of fluid particles by 
turbulence has been first addressed by L.F. Richardson 
[l[ in a pioneering paper in 1926 (see @, \^ for recent 
reviews). The main question is simple and fundamen- 
tal: given a pair of particles released, at time to, at a 
small separation r$ (smaller of the Kolmogorov dissipa- 
tive scale, rf) what is the probability to find the pair at 
any given distance, r, at a later time i? In the case of 
isotropic and homogeneous turbulence (HIT) the proba- 
bility density function (PDF), P(r, t\ro, to), of pair sepa- 
rations depends on the amplitude of r{t) only. Moreover, 
asymptotically it should become independent of the ini- 
tial condition, P(r,t). 

The knowledge of this distribution is of utmost impor- 
tance for many applied studies in geophysics 0-0 , being 
connected also to risk-evaluations in environmental disas- 
ters. Moreover, it constitutes a highly non trivial statis- 
tical problem, due to its intrinsic non-stationary nature 
which connects Lagrangian properties, i.e. velocity fluc- 
tuations evaluated at the particle positions, along their 
whole past history Q. Richardson proposed to model 
the particle separation, at inertial subrange distances 
T] <C r <C Lq, as a diffusive process characterized by an ef- 
fective turbulent diffusivity, Dm c {r) = h ^ , estimated 
empirically to follow a 4/3- law: 



2 dt 



D Ric (r) ~ /3r 4 / 3 



(1) 



Here Lq is the large-scale of the flow and /3 = fcoe 1 / 3 , 
where ko is a dimcnsionless constant and e is the turbu- 
lent kinetic energy dissipation. It is easy Richardson's 
work with the Kolmogorov 1941 theory |8[ by following 
the dimensional estimate 0, H[ : 



D Ric (r) ~ T {r)((6 r v) 2 ) 



(2) 



where r(r) ~ e V 3 ?- 2 / 3 is the eddy-turn-over-timc at 
scale r and {(S r v) 2 } = Coe 2 / 3 r 2 / 3 is the second order 



Eulerian longitudinal structure function. The resulting 
long time growth of the mean squared separation is: 



(r 2 (t))=get 3 



(3) 



g being the Richardson constant uniquely determined in 
terms of fco (lol - [l3l. Sin ce Ej , many experimental and nu- 
merical studies 12J, ll4H16| have focused on the subject, 
inclu ding also extensions to the case of particles with in- 
ertia 17-1^]. There are not doubts that Richardson's 
picture captures some important features of turbulent 
dispersion, in particular concerning events with a typi- 
cal separation of the order of the mean. On the other 
hand, fundamental questions exist about the capability 
of the theory to correctly predict also extremal events, 
i.e. those pairs with separation much larger/smaller than 



1/2 



Such events arc very difficult to measure, due 



to their low probability in absolute terms. Richardson's 
approach can be reinterpreted as the evolution of trac- 
ers in a stochastic Gaussian, homogeneous and isotropic 
velocity field, delta-correlated in time, with two-point 
correlation |2fl|: (5 r Vi(t)5 r Vj(tf)) = 5{t - t')[D\\(r)rirj + 
D±(r)(6ij —fifj)] , where the longitudinal and transverse 
second order structure functions D 11 (r) and D±(r) are 
such that D±(r) = Di\ (r) + l/2rd r D\\ (r), by incompress- 
ibility. Under this assumption, the evolution of P(r, t) is 
closed and local: @,[2o[: 



d t P(r,t) = r^dyD^drPir,!). 



(4) 



Whenever the eddy-diffusivity has a power-law behavior, 
D 11 (r) = Dor* with < £ < 2, the above equation with 
initial condition P(r,to) tx S(r — ro) can be solved in 
terms of special functions [U [22)]. So doing, it admits 
the celebrated asymptotic large time form (independent 
of r ): 

PRic(r,t) oc J 3 exp \-b ( !" 1 ) 1. (5) 
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FIG. 1: Typical time history up to t = 75r v of a bunch emit- 
ted from a source of size ~ rj. Inset: time history for the same 
duration of a bunch emitted in a different location and which 
does not separate. DNS are performed on a cubic fully pe- 
riodic grid at 1024 3 collocation points with a pseudospectral 
code, at a Reynolds- Taylor number, Re\ ~ 300. For further 
details on the numerics see 1121. 



Here b is a constant uniquely determined in terms of Dq 
2l| , which plays the same role of /3 in cq. ([1]). In such ide- 
alized scaling scenario, tracer pairs separate in a explosive 
way forgetting their initial separation : (r 2 (t)) oc t 2 l ( 2 ~£) , 
which reproduces Richardson's expression for £ = 4/3. 
When the turbulent flow is differentiable, i.e. £ = 
2, the PDF takes the log-normal form, P(r, t\ro,0) oc 
cxp{-(log(r/r ) - Xt) 2 /(2A £)}, where A is the first 
Lyapunov exponent and A is connected to fluctuations 
of the strain matrix Q. In the latter case, particles sep- 
arate exponentially and the memory of the initial sepa- 
ration, ro, remains at all times. The rate of separation is 
strongly fluctuating from point-to-point and from time- 
to-time being connected to the fluctuations of the Lya- 
punov exponents @, [23| . 

The behavior of particle pairs in real flows can deviate 
from Richardson distribution because of different rea- 
sons. The four most important are: (i) temporal cor- 
relations of the underlying velocity fluid [H, Ho, HH ; (ii) 
non-Gaussian velocity fluctuations; (iii) ultraviolet (UV) 
effects induced by the dissipation subrange, and (iv) 
infrared (IR) effects induced by the large-scale cut-off. 
These last two are of course connected to finite Reynolds 
effects [25|. 

The goal of this Letter is twofold. First, we want to 
understand and quantify the rate of occurrence of rare 
extreme events, i.e., of pairs that separate much more or 
much less than (r 2 (t)) 1 / 2 . Second, we aim to assess the 



importance of temporal correlations for the pair statistics 
both in general and for extreme events, in particular. 
We have performed a series of Direct Numerical Simu- 
lations (DNS) of HIT seeding the fluid with bunches of 
tracers, emitted in different locations, to reduce spatial 



correlations. Each bunch is emitted within a small region 
of size ~ rj, in puffs of 2 x 10 3 tracer particles each. In a 
single run, there are 256 of such sources, emitting about 
200 puffs with a frequency comparable with the inverse 
of the Kolmogorov time. We performed 10 different runs, 
following a total of 4 x 10 11 pairs, reaching an unprece- 
dented statistics. In Fig. (Q]) we illustrate the complexity 
of the problem. We first notice the abrupt transition in 
the particle dispersion occurring at a time about ~ 10 
after the emission when most of the pairs reaches a rela- 
tive distance of the order of ~ 10rj, and the beginning of 
an explosive separation a la Richardson is observed. At 
at later time, there are many pairs with relative separa- 
tions of the order of the box size ~ IOOO77, even though 
the mean separation is much smaller at those time lags. 
On the contrary, in the inset we show an example of a 
bunch with an anomalous history, due to tracers that 
travel close -at mutual distance of the order of rj -, for 
very long times. This happens when pairs are injected in 
a space location where the underlying fluid has a small 
local stretching rate. Nevertheless, once the pairs reach 
inertial subrange separations, the bunch rapidly expands 
forgetting its initial delay and recovers at large times a 
spread distribution (not shown here). 
To quantify this phenomenology, we show in Fig. (J2J the 
right and left tails of P(r,t) at different time lags, aver- 
aged over all emissions and over all point sources. The 
top panel shows the fastest separation events, with a clear 
exponential-like tail plus a sharp drop at a cut-off separa- 
tion, r c (t), that evolves in time. The cut-off scale r c (t) is 
connected to pairs which are able to separate "very fast" , 
and it is a clear indication of the change in the physics 
governing large excursions. It indicates the presence of a 
finite maximal propagation velocity: it is the signature of 
tracer pairs experiencing a persistently high relative ve- 
locity, which must be limited by the root-mean-squared 
velocity, v rm s [2fi | . To support this statement we show in 
the inset the evolution of r c (t) which is in good agreement 
with a linear behavior obtained using v rms as traveling 
speed. Events beyond the cut-off scale are rare, and can 
be detected with high-statistics only. 
The opposite limit of "very slow" events is also remark- 
able (see bottom panel of Fig. Here we observe a 
bi-modal shape for P(r, t) at almost all times: the left 
tail of pairs with mutual distance r < 77 remains popu- 
lated for time up to ~ 60 — 7QT n , which is of the order 
of the large-scale eddy-turn-over-time, Tl- Pairs emit- 
ted in regions with a small stretching rate tend to stay 
together, reaching only very late the inertial subrange 
separations, and thus never experiencing a Richardson- 
like dispersion. These pairs are governed by a log-normal 
distribution. This is a non-trivial result in tracers disper- 
sion, that can not be brought back to small-scale cluster- 
ing effects as those observed in the dynamics of inertial 
particles [13, HI]. 
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FIG. 2: Top: Log-lin plot of P(r, t) at different times after 
the emission. For selected values of r, we show error bars, 
estimated from the statistical spread of different runs. Inset: 
evolution of the cut-off scale r c (t), the continuous line rep- 
resents the ballistic motion oc v rms t. Bottom: log-log plot 
of P(r,t) for t = (10,20,30,40,50,60,70,90, 120)t„. Black 
squares indicates the peak observed for small separations. 



In Fig. <|3j> we plot the same data of Fig. ([2]) but rcscaled 
in terms of the variable r n (t) = r(t) / (r 2 ^)} 1 / 2 , and com- 
pared against the asymptotic prediction ([SJ. Here, the 
deviations from Richardon's prediction becomes clear, 
showing evident discrepancies at large scales for all times. 
A more stringent test is obtained by showing these same 
PDFs but restricted to the scales in the inertial subrange, 
30r; < r < 300?7 (inset). Clearly ((SJ is not satisfied. Pre- 
vious studies could access events up to r/(r 2 (i)) 1//2 < 3 
only (see Q). Our study improves of four- five order of 
magnitudes (in probability) the intensities of detectable 
events, thus allowing to highlight strong deviations from 
Richardon's shape. Large discrepancies can be measured 
also on the left tails of P(r,t), associated to very slow 
separating pairs (see also below). 

Such departures from the ideal self-similar Richardson 
distribution needs to be better quantified, either in terms 
of finite Reynolds effects (break-up of self-similarity of 
the turbulent eddy diffusivity) or in terms of the ne- 
glected temporal correlations, or both. 
To assess the importance of the former, we have inte- 




FIG. 3: Log-lin plot of P(r n ,t) versus the rescaled variable r„ 
(see text) for t = (20, 30, 40, 60, 90, 120) t>,. The distribution 
P(r n , t) has been divided by a factor r„ to highlight the large 
separation range. Richardson prediction © becomes time- 
independent if rescaled in this way (solid curve). Inset: same 
PDFs plotted only for separations r n that, at any time lag, 
belong to the inertial subrange. 



grated the Richardson equation ^ using an effective 
eddy-diffusivity D e * f (r) which improves ^ by including 
UV and IR cut-offs. Depending whether the separation 
distance falls in the viscous, inertial or integral range of 
scales, we must then have : 



Df f (r) 



D?/ f (r) 



r r <C ?y 

r 4/3 < r <C L 

const, r > Iq- 



(6) 



A widely used fitting formula that reproduces well the 
Eulerian data, and that matches the expected UV and 
IR scaling for both r(r) and ((S r v) 2 }, is obtained by fol- 
lowing [29( : 



((Srvf 



CO 



((r/77) 2 +Ci)( 2 / 3 ) 



1 -1/3 



l + c 2 ( — f 

-^0 



(7) 



supplemented with a similar expression for the eddy turn 
overtime: r(r) = {{r/7l)2 + Ci) - 1/3 [l + d 2 (r/L ) 2 ] ~ 2/3 . 
The dimcnsionless parameters Cq , C\ , c 2 are extracted 
from the Eulerian statistics, while the parameter d 2 is 
fixed such as to correctly reproduce the evolution of 
the mean square separation, (r 2 (t)), over a time range 
t v < t < Tl, see Fig. |4] Note that the hypothesis 
of Gaussian statistics still (implicitly) holds in this ap- 
proach, since the velocity field distribution is fixed in 
terms of the second order moment only. 

The solution, P e ff(r,t), to the diffusive equation (|4]) 
using D e ff(r) is shown in Fig. [5] and compared with 
the DNS data. Despite the excellent agreement between 
{r 2 (t)) obtained from DNS and the one obtained using 
the stochastic model ^ and ©, the far tails are still 
different. Self-similarity is broken by the introduction of 
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FIG. 4: Log-log plot of (r 2 (t)) from DNS data, and from 
the diffusive evolution with eddy diffusivity ([6]). Inset: time 
evolution of the relative probability to observe a large excur- 
sion, 7>(A,t) (right y-scale), or small excursion I<(A,f) (left 
y-scale) for A = 3. 



UV and IR cutoffs in Eqn. ©, and therefore P e ff(r, t) no 
longer rescales at different times: this is at variance with 
Richardson idealized picture, and clearly goes in the di- 
rection of real turbulent flows. For large separations the 
agreement with the DNS data is qualitatively better but 
still quantitatively off, particularly when focusing on the 
sharp change at the cut-off scale r c (i) which is still ab- 
sent in the evolution given by (j6)). This is a key-point 
showing that to reproduce the observed drop in the PDF 
for large excursions, it is not enough to impose a sat- 
uration of eddy-diffusivity for large r. The behavior of 
pair dispersion must then be crucially dependent on the 
nature of temporal correlations of the Lagrangian tur- 
bulent velocity in the inertial subrange. Similarly, if we 
focus on the small separation tail, P e f / (r, t) presents a 
slowly evolving peak for r < r\ but quantitative agree- 
ment is again not satisfactory yet (inset of Fig. [5]). Wc 
understand such departure as the effects of assuming a 
Gaussian velocity statistics underlying the evolution of 
P e f f (r, t) , which is blatantly wrong because of turbulent 
small-scale intermittency. To further quantify the depar- 
ture of the modified Richardson description from the real 
data we measured the cumulative probability to have a 
couple at large separation r*(A) = A(r 2 (t)) 1 / 2 normalized 
with the same quantity evaluated from Eqn. ([3]) using 
©, namely /> (A, t) = f™ (x) drP(r,t)/ /~ drP eff (r,t); 
similarly, to evaluate the differences for small separa- 
tion events, we use r*(A) = 1/A(r 2 (i)} 1//2 and define: 

I<(X,t) = jf W drP(r,t)/ J r * drP eff {r,t). The results 
are shown in the inset of Fig. [4] for A = 3. Concern- 
ing 7<(A,i), it shows that the evolution using P e ff(r,t) 
underestimate of a factor 2 at small time-lags, ~ 20^ 
the importance of the small-scale trapping properties, i.e. 
does not capture the strong intermittency of the regions 
where we have a small stretching rate. Only for very large 




FIG. 5: Lin-log plot of P e //(r, t) as obtained from the inte- 
gration of © with D eff (r) (dashed) and the DNS data (solid 
line). Inset dog- log plot to highlight the slowest events (left 
tail). 

times ~ IOOt,, the left tail becomes comparable with the 
real ones. Concerning 7> (A, t) we measure first an un- 
derestimate of large separation events and later a strong 
overestimate, i.e. the delta-correlated approach behind 
the evolution of P e ff(r,t) does not capture the presence 
of the cut-off scale r c (t). 

To directly quantify the importance of temporal corre- 
lations along tracer trajectories, we compare in Fig. ((5]) 
Lagrangian longitudinal velocity increments conditioned 
on particle separation distance r : 

S L ag(r,t) = ((<J r (t)Vi ' h)\r(t) = r) , (8) 

with its Eulerian equivalent : 

SEui(r) = (\S r Vi ■ fi\) . (9) 

The figure shows that pairs that have separated fast have 
also a typical velocity difference higher than the Eulerian 
one, measured on the whole volume without condition- 
ing on the history of particles. This is in agreement with 
what found before: intense separation events feel highly 
persistent velocity increments along their time history, 
and they cannot be trivially estimated starting from Eu- 
lerian statistics. Figure ^ also shows the opposite trend: 
pairs that have not separated much, at any given time, 
have Lagrangian velocity increments smaller than the Eu- 
lerian ones. 

In conclusion, we have presented the first numerical 
study at very high statistics meant to explore rare events 
in turbulent dispersion. Both "fast" and "slow" separa- 
tions show an important departure from the Richardson- 
like inertial and idealized behavior of Eqn. ([5]). A step 
forward toward a better modelization is obtained by 
maintaining the assumption of time delta-correlation and 
Gaussian statistics but introducing an improved effec- 
tive eddy-diffusivity kernel, which takes correctly into 
account both the viscous and integral scale physics. Do- 
ing so, a better qualitative agreement with real data is 
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FIG. 6: Comparison between conditional Lagrangian (thin 
lines) and Eulerian (thick line) velocity increment moment. 



observed. To further progress we will need to relax the 
Gaussianity assumption for the small-separation and the 
delta-correlated assumption for separations of the order 
of the integral scale. An attempt following the latter di- 
rection by using a Langevin-process for the evolution of 
the relative particles velocity, with different correlation 
times at different scales, has been recently proposed by 
[30j . Alternatively, a memory kernel can be introduced 
in order to overcome the presence of unphysical infinite 
speed events 24j . characterizing pairs that move in a 
delta-correlated velocity field. For example, Id model 



with telegraph equations shows similar maximum speed 
events [26[. 
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